dweibull (Double Weibull) Distribution#
This notebook is a self-contained, math-first tour of SciPy’s scipy.stats.dweibull distribution.
Goals
Understand the shape parameter and how it controls peakedness, bimodality, and tail behavior.
Derive key results (PDF/CDF, moments, likelihood) with clean substitutions.
Implement NumPy-only sampling and validate it against SciPy.
See how
dweibullshows up in testing, Bayesian inference, and generative/noise models.
Throughout, we use the standardized distribution (default loc=0, scale=1) unless stated otherwise.
import numpy as np
import plotly.express as px
import plotly.graph_objects as go
import os
import plotly.io as pio
from scipy.special import gamma
from scipy.stats import dweibull, kstest
pio.templates.default = "plotly_white"
pio.renderers.default = os.environ.get("PLOTLY_RENDERER", "notebook")
np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(42)
1) Title & Classification#
Name (SciPy):
dweibull(double Weibull distribution)Type: Continuous
Support: \(x \in (-\infty, \infty)\)
Parameter space (standard form): shape \(c > 0\)
SciPy parameterization:
dweibull(c, loc=0, scale=1)with\(c > 0\) (shape)
\(\text{loc} \in \mathbb{R}\) (location shift)
\(\text{scale} > 0\) (scale stretch)
2) Intuition & Motivation#
A simple generative story#
A convenient way to think about dweibull is:
Draw a magnitude \(Y \ge 0\) from a (one-sided) Weibull distribution.
Flip a fair coin for the sign \(S \in \{-1, +1\}\).
Set \(X = S\,Y\).
This immediately explains why the distribution is symmetric about 0 and why \(|X|\) is Weibull.
What it models#
dweibull is useful when you want a symmetric distribution whose shape can morph between:
sharp peak at 0 with heavy (stretched-exponential) tails (\(0 < c < 1\))
Laplace / double-exponential (\(c = 1\))
bimodal shapes with a dip at 0 (\(c > 1\))
The key surprise is the last bullet: for \(c>1\), the PDF at 0 becomes zero, and the distribution has two symmetric modes away from 0.
Real-world use cases#
Signed magnitudes: deviations that come with a size (Weibull-like) and a random sign (e.g., anomaly sizes, symmetric measurement deviations).
Flexible error/noise models: use \(c\) to tune tail heaviness vs concentration near 0 (especially for \(c\le 1\)).
Bimodal symmetric data: when values tend to avoid 0 but cluster around \(\pm m\) for some magnitude.
Relations to other distributions#
If \(X \sim \texttt{dweibull}(c)\), then \(|X|\) is Weibull with the same shape parameter \(c\).
\(c=1\) gives the Laplace distribution: \(f(x)=\tfrac12 e^{-|x|}\).
\(c=2\) implies \(|X|\) is Rayleigh (with a particular scale), so
dweibullbecomes a symmetric “signed Rayleigh magnitude” model.
3) Formal Definition#
PDF (standardized)#
For shape \(c>0\), the probability density function is
CDF (standardized)#
The CDF has a clean piecewise form. For \(x<0\),
and for \(x\ge 0\),
Location/scale form#
SciPy’s loc and scale apply the standard transformation
Then
Quantile function (PPF)#
Because the CDF is explicit, inverse-CDF sampling is easy. For \(q\in(0,1)\),
def _validate_params(c: float, scale: float) -> None:
if not np.isfinite(c) or c <= 0:
raise ValueError(f"shape c must be > 0, got {c!r}")
if not np.isfinite(scale) or scale <= 0:
raise ValueError(f"scale must be > 0, got {scale!r}")
def dweibull_pdf(x, c: float, loc: float = 0.0, scale: float = 1.0):
# NumPy implementation of the PDF (with loc/scale).
_validate_params(c, scale)
x = np.asarray(x, dtype=float)
z = (x - loc) / scale
az = np.abs(z)
return (c / (2 * scale)) * np.power(az, c - 1) * np.exp(-np.power(az, c))
def dweibull_logpdf(x, c: float, loc: float = 0.0, scale: float = 1.0):
# NumPy implementation of log-PDF (handles x==loc explicitly).
_validate_params(c, scale)
x = np.asarray(x, dtype=float)
z = (x - loc) / scale
az = np.abs(z)
out = np.empty_like(az)
zero = az == 0
# At z=0: pdf(0) is 0 if c>1, 1/(2*scale) if c=1, and +inf if c<1.
if c > 1:
out[zero] = -np.inf
elif np.isclose(c, 1.0):
out[zero] = -np.log(2 * scale)
else:
out[zero] = np.inf
nz = ~zero
out[nz] = (
np.log(c)
- np.log(2 * scale)
+ (c - 1) * np.log(az[nz])
- np.power(az[nz], c)
)
return out
def dweibull_cdf(x, c: float, loc: float = 0.0, scale: float = 1.0):
# NumPy implementation of the CDF using expm1 for precision near 0.
_validate_params(c, scale)
x = np.asarray(x, dtype=float)
z = (x - loc) / scale
out = np.empty_like(z)
neg = z < 0
pos = ~neg
zneg = -z[neg]
zpos = z[pos]
out[neg] = 0.5 * (1.0 + np.expm1(-np.power(zneg, c)))
out[pos] = 0.5 - 0.5 * np.expm1(-np.power(zpos, c))
return out
def dweibull_ppf(q, c: float, loc: float = 0.0, scale: float = 1.0):
# Quantile function (inverse CDF).
_validate_params(c, scale)
q = np.asarray(q, dtype=float)
out = np.empty_like(q)
out[q <= 0] = -np.inf
out[q >= 1] = np.inf
mid = (q > 0) & (q < 1)
qmid = q[mid]
left = qmid < 0.5
right = ~left
out_mid = np.empty_like(qmid)
out_mid[left] = -np.power(-np.log(2.0 * qmid[left]), 1.0 / c)
out_mid[right] = np.power(-np.log(2.0 * (1.0 - qmid[right])), 1.0 / c)
out[mid] = loc + scale * out_mid
return out
# Quick sanity check: CDF(PPF(q)) ≈ q
c_test = 1.3
q_grid = np.linspace(1e-6, 1 - 1e-6, 2000)
x_from_q = dweibull_ppf(q_grid, c=c_test)
q_back = dweibull_cdf(x_from_q, c=c_test)
max_err = np.max(np.abs(q_back - q_grid))
max_err
1.1102230246251565e-16
4) Moments & Properties#
Absolute moments (key identity)#
For the standardized distribution \(Z\sim\texttt{dweibull}(c)\), a very useful identity is
In particular, all positive moments exist for any \(c>0\).
Mean, variance, skewness, kurtosis#
Because the PDF is symmetric, all odd central moments are 0 (when they exist). For \(c>0\):
Mean: \(\mathbb{E}[Z]=0\) and \(\mathbb{E}[X]=\text{loc}\).
Variance: \(\operatorname{Var}(Z)=\Gamma\!\left(1+\frac{2}{c}\right)\), so \(\operatorname{Var}(X)=\text{scale}^2\,\Gamma\!\left(1+\frac{2}{c}\right)\).
Skewness: 0.
Excess kurtosis:
MGF / characteristic function#
The characteristic function \(\varphi(t)=\mathbb{E}[e^{itZ}]\) always exists. Because of symmetry, \(\varphi(t)=\mathbb{E}[\cos(tZ)]\).
The MGF \(M(t)=\mathbb{E}[e^{tZ}]\) depends on \(c\):
\(c>1\): exists for all real \(t\) (tails decay faster than exponential).
\(c=1\): exists only for \(|t|<1\) (Laplace case).
\(0<c<1\): diverges for any \(t\ne 0\) (stretched-exponential tails).
A useful analytic representation is the even-moment series (when it converges):
Entropy#
The differential entropy for the standardized distribution is
where \(\gamma\approx 0.57721\) is the Euler–Mascheroni constant. With scaling, \(H(X)=H(Z)+\ln(\text{scale})\).
Modes#
For \(c\le 1\), the distribution is unimodal with a peak at 0 (in fact, the PDF is infinite at 0 when \(c<1\)). For \(c>1\), the PDF at 0 is 0 and there are two modes at
def dweibull_theoretical_stats(c: float, loc: float = 0.0, scale: float = 1.0):
# Return mean, var, skewness, excess kurtosis for dweibull(c, loc, scale).
_validate_params(c, scale)
mean = loc
var = (scale**2) * gamma(1.0 + 2.0 / c)
# Symmetry => skewness = 0
skew = 0.0
m2 = gamma(1.0 + 2.0 / c)
m4 = gamma(1.0 + 4.0 / c)
excess_kurtosis = m4 / (m2**2) - 3.0
return mean, var, skew, excess_kurtosis
def dweibull_entropy(c: float, scale: float = 1.0):
# Differential entropy for dweibull(c, loc=0, scale).
_validate_params(c, scale)
return 1.0 - np.log(c) + np.log(2.0 * scale) + np.euler_gamma * (1.0 - 1.0 / c)
c_demo = 0.7
mean_th, var_th, skew_th, kurt_th = dweibull_theoretical_stats(c_demo)
H_th = dweibull_entropy(c_demo)
mean_sp, var_sp, skew_sp, kurt_sp = dweibull.stats(c_demo, moments="mvsk")
H_sp = dweibull.entropy(c_demo)
(
np.array([mean_th, var_th, skew_th, kurt_th]),
np.array([mean_sp, var_sp, skew_sp, kurt_sp]),
float(H_th),
float(H_sp),
)
(array([ 0. , 5.029145, 0. , 13.777367]),
array([ 0. , 5.029145, 0. , 13.777367]),
1.802443982398021,
1.802443982398021)
5) Parameter Interpretation (How Shape Changes)#
Shape parameter \(c\)#
The single shape parameter controls multiple behaviors at once:
Near 0: the factor \(|x|^{c-1}\) decides what happens at the origin.
\(c<1\): \(|x|^{c-1}\to\infty\) → infinite spike at 0.
\(c=1\): finite value at 0 (Laplace).
\(c>1\): \(|x|^{c-1}\to 0\) → density drops to 0 at 0 (bimodal).
Tails: \(\exp(-|x|^c)\) controls tail decay.
Smaller \(c\) → heavier (slower) tail decay.
Larger \(c\) → lighter (faster) tail decay.
loc and scale#
locshifts the distribution left/right (median and mean move toloc).scalestretches the distribution; variance scales like \(\text{scale}^2\).
# PDF shapes for different c
x = np.linspace(-4, 4, 2000)
c_values = [0.5, 0.8, 1.0, 1.5, 3.0]
fig = go.Figure()
for c in c_values:
y = dweibull_pdf(x, c)
# For c<1, the PDF spikes to +inf at 0; clip just for plotting.
finite = np.isfinite(y)
if np.any(finite):
cap = np.nanquantile(y[finite], 0.995)
y_plot = np.clip(y, 0, cap)
else:
y_plot = y
fig.add_trace(go.Scatter(x=x, y=y_plot, mode="lines", name=f"c={c}"))
fig.update_layout(
title="dweibull PDF shapes (clipped near the spike for c<1)",
xaxis_title="x",
yaxis_title="pdf(x)",
)
fig.show()
# CDF shapes for the same c values
x = np.linspace(-4, 4, 2000)
fig = go.Figure()
for c in c_values:
fig.add_trace(go.Scatter(x=x, y=dweibull_cdf(x, c), mode="lines", name=f"c={c}"))
fig.update_layout(
title="dweibull CDF shapes",
xaxis_title="x",
yaxis_title="cdf(x)",
)
fig.show()
6) Derivations#
(a) Expectation#
For the standardized distribution \(Z\) the PDF is symmetric: \(f(z)=f(-z)\). Provided \(\mathbb{E}[|Z|]<\infty\) (true for all \(c>0\)),
With location, \(X=\text{loc}+\text{scale}Z\), we get \(\mathbb{E}[X]=\text{loc}\).
(b) Variance via the Gamma function#
Compute the absolute moment for \(k> -c\):
Substitute \(u=z^c\Rightarrow z=u^{1/c}\) and \(dz=\tfrac1c u^{1/c-1}du\):
For the variance, take \(k=2\):
(c) Likelihood and log-likelihood#
Assume i.i.d. observations \(x_1,\dots,x_n\) from the standardized model (\(\text{loc}=0\), \(\text{scale}=1\)). The likelihood for \(c\) is
The log-likelihood is
Differentiating gives a score equation (no closed-form MLE in general):
This is typically solved numerically (as SciPy does under the hood).
# Example: visualize the log-likelihood over c (standardized case)
x_data = dweibull.rvs(0.8, size=4000, random_state=rng)
def loglike_c(c: float) -> float:
if c <= 0:
return -np.inf
return float(np.sum(dweibull_logpdf(x_data, c)))
c_grid = np.linspace(0.2, 4.0, 200)
ll = np.array([loglike_c(c) for c in c_grid])
fig = px.line(x=c_grid, y=ll, labels={"x": "c", "y": "log-likelihood"}, title="Log-likelihood vs c")
fig.show()
c_hat = float(c_grid[np.argmax(ll)])
c_hat
0.7919597989949749
7) Sampling & Simulation (NumPy-only)#
Inverse transform sampling#
Using the PPF, we can sample with a single uniform random variable, but an even simpler implementation uses the sign + magnitude story:
Sample \(U\sim\text{Unif}(0,1)\) and set \(Y = (-\ln U)^{1/c}\). (This is Weibull sampling.)
Sample an independent sign \(S\in\{-1,+1\}\) with \(\mathbb{P}(S=1)=1/2\).
Return \(X = \text{loc} + \text{scale}\cdot S Y\).
This is NumPy-only and avoids SciPy entirely.
def dweibull_rvs_numpy(
c: float,
loc: float = 0.0,
scale: float = 1.0,
size: int | tuple[int, ...] = 1,
rng: np.random.Generator | None = None,
):
# Draw random samples from dweibull using NumPy only.
_validate_params(c, scale)
if rng is None:
rng = np.random.default_rng()
u = rng.random(size)
u = np.clip(u, np.finfo(float).tiny, 1.0)
# Magnitude ~ Weibull(shape=c, scale=1): Y = (-log U)^(1/c)
y = np.power(-np.log(u), 1.0 / c)
# Random sign
s = np.where(rng.random(size) < 0.5, -1.0, 1.0)
return loc + scale * s * y
# Quick validation: sample moments vs theory
c_sim = 0.8
n = 100_000
samples = dweibull_rvs_numpy(c_sim, size=n, rng=rng)
mean_mc = float(np.mean(samples))
var_mc = float(np.var(samples))
mean_th, var_th, *_ = dweibull_theoretical_stats(c_sim)
(mean_mc, var_mc, float(mean_th), float(var_th))
(-0.014103101477437954, 3.3469432497508005, 0.0, 3.323350970447843)
8) Visualization#
We’ll visualize three things for a chosen parameter set:
the PDF
the CDF
Monte Carlo samples compared to the theoretical PDF
c_vis = 0.8
loc_vis = 0.0
scale_vis = 1.2
x = np.linspace(-5, 5, 3000)
pdf_np = dweibull_pdf(x, c_vis, loc=loc_vis, scale=scale_vis)
cdf_np = dweibull_cdf(x, c_vis, loc=loc_vis, scale=scale_vis)
pdf_sp = dweibull.pdf(x, c_vis, loc=loc_vis, scale=scale_vis)
cdf_sp = dweibull.cdf(x, c_vis, loc=loc_vis, scale=scale_vis)
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=pdf_np, mode="lines", name="NumPy pdf"))
fig.add_trace(go.Scatter(x=x, y=pdf_sp, mode="lines", name="SciPy pdf", line=dict(dash="dash")))
fig.update_layout(title="PDF: NumPy vs SciPy", xaxis_title="x", yaxis_title="pdf(x)")
fig.show()
fig = go.Figure()
fig.add_trace(go.Scatter(x=x, y=cdf_np, mode="lines", name="NumPy cdf"))
fig.add_trace(go.Scatter(x=x, y=cdf_sp, mode="lines", name="SciPy cdf", line=dict(dash="dash")))
fig.update_layout(title="CDF: NumPy vs SciPy", xaxis_title="x", yaxis_title="cdf(x)")
fig.show()
# Monte Carlo histogram vs theoretical PDF
n = 20_000
x_samp = dweibull_rvs_numpy(c_vis, loc=loc_vis, scale=scale_vis, size=n, rng=rng)
hist = px.histogram(
x=x_samp,
nbins=80,
histnorm="probability density",
opacity=0.5,
title="Samples (histogram) vs theoretical PDF",
labels={"x": "x"},
)
pdf_line = go.Scatter(x=x, y=pdf_sp, mode="lines", name="theoretical pdf")
fig = go.Figure(hist.data)
fig.add_trace(pdf_line)
fig.update_layout(xaxis_title="x", yaxis_title="density")
fig.show()
9) SciPy Integration (scipy.stats.dweibull)#
SciPy offers a full suite of distribution methods:
pdf,logpdf,cdf,ppf,rvsstatsfor momentsentropyfitfor parameter estimation (numerical)
We’ll generate synthetic data, fit the parameters, and overlay the fitted PDF.
# Synthetic data
c_true, loc_true, scale_true = 0.9, -0.3, 1.4
x_obs = dweibull.rvs(c_true, loc=loc_true, scale=scale_true, size=5_000, random_state=rng)
# Fit all parameters
c_fit, loc_fit, scale_fit = dweibull.fit(x_obs)
(c_fit, loc_fit, scale_fit)
(0.8918390627792722, -0.2964891361570203, 1.3802126275005553)
x_grid = np.linspace(np.quantile(x_obs, 0.001), np.quantile(x_obs, 0.999), 2000)
pdf_true = dweibull.pdf(x_grid, c_true, loc=loc_true, scale=scale_true)
pdf_fit = dweibull.pdf(x_grid, c_fit, loc=loc_fit, scale=scale_fit)
hist = px.histogram(
x=x_obs,
nbins=80,
histnorm="probability density",
opacity=0.4,
title="SciPy fit: true vs fitted PDF",
labels={"x": "x"},
)
fig = go.Figure(hist.data)
fig.add_trace(go.Scatter(x=x_grid, y=pdf_true, mode="lines", name="true pdf"))
fig.add_trace(go.Scatter(x=x_grid, y=pdf_fit, mode="lines", name="fitted pdf", line=dict(dash="dash")))
fig.update_layout(xaxis_title="x", yaxis_title="density")
fig.show()
10) Statistical Use Cases#
(a) Hypothesis testing (goodness-of-fit)#
A common workflow is:
fit parameters
test whether the fitted distribution plausibly generated the data
A classic tool is the Kolmogorov–Smirnov (KS) test.
Caution: if you fit parameters on the same data you test, the KS p-value is only approximate (the null distribution changes). Still, it’s a useful diagnostic.
(b) Bayesian modeling#
The log-likelihood \(\ell(c)\) makes it easy to do Bayesian inference for \(c\) with a prior (e.g., Gamma prior). We’ll do a simple grid posterior example with known loc=0, scale=1.
(c) Generative modeling#
You can use dweibull as a drop-in noise distribution (especially for \(c\le 1\)) to generate data with heavier tails than a Gaussian and a different near-zero behavior.
# (a) KS test using fitted parameters (approximate when parameters are estimated)
D, p_value = kstest(x_obs, "dweibull", args=(c_fit, loc_fit, scale_fit))
(D, p_value)
(0.015158632683738738, 0.19873182866725814)
# (b) Simple Bayesian inference for c (assuming loc=0, scale=1 known)
# Generate standardized data
c_true_bayes = 0.75
x_bayes = dweibull.rvs(c_true_bayes, size=1500, random_state=rng)
c_grid = np.linspace(0.2, 4.0, 600)
# Gamma prior on c: shape α, rate β
alpha, beta = 2.0, 1.0
log_prior = (alpha - 1) * np.log(c_grid) - beta * c_grid # constants omitted
log_like = np.array([np.sum(dweibull_logpdf(x_bayes, c)) for c in c_grid])
log_post = log_like + log_prior
log_post -= np.max(log_post)
post_unnorm = np.exp(log_post)
# Normalize
Z = np.trapz(post_unnorm, c_grid)
post = post_unnorm / Z
# Posterior mean and MAP
c_map = float(c_grid[np.argmax(post)])
c_mean = float(np.trapz(c_grid * post, c_grid))
# 95% credible interval via numerical CDF
cdf = np.cumsum((post[:-1] + post[1:]) / 2 * np.diff(c_grid))
cdf = np.concatenate([[0.0], cdf])
c_lo = float(np.interp(0.025, cdf, c_grid))
c_hi = float(np.interp(0.975, cdf, c_grid))
(c_true_bayes, c_map, c_mean, (c_lo, c_hi))
(0.75,
0.7455759599332219,
0.7472222579871374,
(0.7182876018799911, 0.7762668748027535))
fig = go.Figure()
fig.add_trace(go.Scatter(x=c_grid, y=post, mode="lines", name="posterior"))
fig.add_vline(x=c_true_bayes, line_dash="dash", line_color="black", annotation_text="true c")
fig.add_vline(x=c_map, line_dash="dot", line_color="royalblue", annotation_text="MAP")
fig.update_layout(title="Posterior over c (loc=0, scale=1 assumed)", xaxis_title="c", yaxis_title="density")
fig.show()
# (c) Generative modeling example: a smooth signal + different noise models
t = np.linspace(0, 1, 400)
y_true = np.sin(2 * np.pi * t)
sigma = 0.25
noise_gauss = rng.normal(0.0, sigma, size=t.size)
noise_dw = dweibull_rvs_numpy(c=0.7, scale=sigma, size=t.size, rng=rng)
y_gauss = y_true + noise_gauss
y_dw = y_true + noise_dw
fig = go.Figure()
fig.add_trace(go.Scatter(x=t, y=y_true, mode="lines", name="true signal", line=dict(color="black")))
fig.add_trace(go.Scatter(x=t, y=y_gauss, mode="markers", name="Gaussian noise", opacity=0.6))
fig.add_trace(go.Scatter(x=t, y=y_dw, mode="markers", name="dweibull noise (c=0.7)", opacity=0.6))
fig.update_layout(title="Same signal, different noise distributions", xaxis_title="t", yaxis_title="y")
fig.show()
# Compare residual distributions
residuals = {
"Gaussian": noise_gauss,
"dweibull (c=0.7)": noise_dw,
}
fig = px.histogram(
x=np.concatenate(list(residuals.values())),
color=np.repeat(list(residuals.keys()), repeats=[t.size, t.size]),
nbins=70,
barmode="overlay",
histnorm="probability density",
opacity=0.5,
title="Residual distributions",
labels={"x": "residual"},
)
fig.show()
11) Pitfalls#
Parameter validity: \(c>0\),
scale>0. Invalid values should error early.Zeros in data: for \(c>1\), the PDF at 0 is exactly 0; if your data contains many exact zeros (rounding/quantization), the likelihood can behave strangely. For \(c<1\), the PDF is infinite at 0, so exact zeros can dominate fits.
Bimodality for \(c>1\): this is often unexpected if you think of the model as “noise around 0”.
MGF nonexistence: for \(c\le 1\), the MGF does not exist for all \(t\) (Laplace has a finite strip; \(c<1\) diverges for any nonzero \(t\)).
Numerical stability: prefer
logpdfin optimization; for large \(|x|^c\) the PDF underflows to 0 (fine), but products of PDFs can underflow without logs.Fitting:
dweibull.fitis numerical and can be sensitive; consider fixinglocif you know the center, or providing good initial guesses in custom optimization.
12) Summary#
dweibullis a continuous, symmetric distribution on \(\mathbb{R}\) with shape parameter \(c>0\).Its PDF \(\propto |x|^{c-1}e^{-|x|^c}\) creates three regimes: spike at 0 (\(c<1\)), Laplace (\(c=1\)), and bimodal (\(c>1\)).
Key identity: \(\mathbb{E}[|Z|^k]=\Gamma(1+k/c)\), giving closed-form variance and kurtosis.
Sampling is simple via inverse transform / sign + Weibull magnitude, and SciPy provides
pdf/cdf/rvs/fitutilities.In practice,
dweibullcan be a flexible tool for diagnostics, Bayesian inference over shape, and generative noise modeling when Gaussian assumptions are not appropriate.
References
SciPy docs:
scipy.stats.dweibull(notes include the defining PDF)Standard Gamma function identities for Weibull moments